## SI figure 2
rm(list=ls())

## Load HadCRUT temperature time series
path <- "~/Documents/Projects/IAM_Stochastic/" # Directory where HadCRUT time series is stored.
HadCRUT <- read.table(paste(path, "HadCRUT.4.5.0.0.annual_ns_avg.txt", sep=""), header=FALSE)
HadCRUT <- HadCRUT[,1:2]
names(HadCRUT) <- c("year","temp")
sf = 2/3 # Smoothing factor. A larger value returns a smoother trend line.
HadCRUT.detrended <- HadCRUT$temp - lowess(HadCRUT$year, HadCRUT$temp, f=sf)$y 

plot(HadCRUT, type="l")
lines(lowess(HadCRUT$year, HadCRUT$temp, f=sf))
lines(HadCRUT$year,lm(HadCRUT$temp ~ HadCRUT$year)$coef[1] + lm(HadCRUT$temp ~ HadCRUT$year)$coef[2]*HadCRUT$year)

## Mean of the HadCRUT time series
HadCRUT.mean <- mean(HadCRUT$temp)

## Slope of the HadCRUT time series
HadCRUT.beta <- lm(HadCRUT$temp ~ HadCRUT$year)$coef[2]

## Variance of the detrended HadCRUT time series
HadCRUT.var <- sum((HadCRUT.detrended)^2)/(nrow(HadCRUT)-1)

## One-lag autocorrelation in the HadCRUT time series
HadCRUT.phi <- cor(HadCRUT$temp[-1], HadCRUT$temp[-nrow(HadCRUT)])
sf = 1 # Smoothing factor. A larger value returns a smoother trend line.
HadCRUT.detrended <- HadCRUT$temp - lowess(HadCRUT$year, HadCRUT$temp, f=sf)$y 
HadCRUT.phi.detrended <- cor(HadCRUT.detrended[-1], HadCRUT.detrended[-nrow(HadCRUT)])

plot(HadCRUT, type="l")
lines(lowess(HadCRUT$year, HadCRUT$temp, f=sf))
lines(HadCRUT$year,lm(HadCRUT$temp ~ HadCRUT$year)$coef[1] + lm(HadCRUT$temp ~ HadCRUT$year)$coef[2]*HadCRUT$year)


## Estimate variance of the Hasselman runs with different (sigmaQ,lambda,Ceff) combinations
files <- list.files(path=path, pattern="Historical_")
d <- data.frame()
for (i in 1:length(files)){
	ds <- read.table(paste(path, files[i], sep=""), sep=",", header=FALSE)
	ds$param <- files[i]
	d <- rbind(d, ds)
}

## Parameter values for sigmaQ and lambda for which the Hasselman model has been run
d$ECS <- as.numeric(sapply(strsplit(d$param, "_"), "[", 2))
d$sigmaQ <- as.numeric(sapply(strsplit(d$param, "_"), "[", 3))
d$Ceff <- as.numeric(sapply(strsplit(d$param, "_"), "[", 4))

## Ensemble size
nn = ncol(d) - 5
names(d) <- c("year",paste("X",seq(1:nn),sep=""), "param","ECS","sigmaQ","Ceff")
d$year <- d$year + 1850

## Length of time series
n=length(unique(d$year))

## Rebase each time series to 1961-1990 average, as HadCRUT, and compute the average variance for each (sigmaQ,lambda) pair.
x <- data.frame()
for (j in unique(d$param)) {
	t <- d[d$param==j,]
	year <- t$year
	ECS <- unique(t$ECS)
	sigmaQ <- unique(t$sigmaQ)
	Ceff <- unique(t$Ceff)

	### Rebase
	t.base <- colMeans(t[t$year %in% 1961:1990,!(names(t) %in% c("year","param","sigmaQ","ECS","Ceff"))])
	t.base <- matrix(rep(t.base,n),ncol=nn,nrow=n,byrow=TRUE)
	t <- t[,!(names(t) %in% c("year","param","sigmaQ","ECS","Ceff"))]
	t <- t - t.base

	### Compute descriptive statistics for Hasselman runs
	means <- mean(as.matrix(t))				# Compute mean
	betas <- rep(NA,nn)
	phis <- rep(NA,nn)
	phis.detrended <- rep(NA,nn)
	for (i in 1:nn) {
		betas[i] <- lm(t[,i] ~ year)$coef[2] # Compute slope
		phis[i] <- cor(t[-1,i], t[-n,i])	 # Compute 1-lag correlation
		t[,i] <- t[,i] - lowess(year, t[,i], f=.2)$y # Detrend
		phis.detrended[i] <- cor(t[-1,i], t[-n,i])	 # Compute detrended 1-lag correlation
	}
	t <- t^2
	t <- colSums(t)
	t <- t/(length(year)-1) # Compute detrended variances
	out <- c(means,mean(t),mean(betas),mean(phis),mean(phis.detrended),ECS,sigmaQ,Ceff)
	x <- rbind(x,out)
}
names(x) <- c("mean","sigmaT","beta","phi","phi.detrended","ECS","sigmaQ","Ceff")

quartz(width = 10, height = 3, type = "pdf", "SI_figure2", file = "SI_figure2.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.5,0,0,.3), # global margins in inches (bottom, left, top, right)
	mai=c(0,.3,.5,.1)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1,2,3), 1, 3, byrow = TRUE), widths=c(3,3,3), heights=c(3,3,3))

# Panel 1
plot(HadCRUT$year,HadCRUT$temp+1,type="l",ylim=c(0,5), bty="L",yaxt="n",xlab="",ylab="",yaxs="i")
axis(side=2,at=seq(-1,10,1),las=2)
Ceff = .8*10^9
ECS = 3
sigmaQ = .625*10^8

lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 4]+1.4, col="grey")
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 27]+2.4, col="grey")
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 30]+3.4, col="grey")
text(bquote("C=" * .(Ceff/10^9) * "x" * 10^9 * "; " * lambda * "=" * .(round(3.7/ECS,2)) * "; " * sigma[Q] * "=" * .(sigmaQ/10^8) * "x" * 10^8), x=1850, y=4.7, pos=4)

# Panel 2
plot(HadCRUT$year,HadCRUT$temp+1,type="l",ylim=c(0,5), bty="L",yaxt="n",xlab="",ylab="",yaxs="i")
axis(side=2,at=seq(-1,10,1),las=2)
Ceff = 1*10^9
ECS = 3
sigmaQ = .9375*10^8
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 98]+1.4, col="grey")
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 71]+2.4, col="grey")
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 51]+3.2, col="grey")
text(bquote("C=" * .(Ceff/10^9) * "x" * 10^9 * "; " * lambda * "=" * .(round(3.7/ECS,2)) * "; " * sigma[Q] * "=" * .(sigmaQ/10^8) * "x" * 10^8), x=1850, y=4.7, pos=4)
text(bquote("C=" * .(Ceff/10^9) * "x" * 10^9 * "; " * lambda * "=" * .(round(3.7/ECS,2)) * "; " * sigma[Q] * "=" * .(sigmaQ/10^8) * "x" * 10^8), x=1850, y=4.7, pos=4)

# Panel 3
plot(HadCRUT$year,HadCRUT$temp+1,type="l",ylim=c(0,5), bty="L",yaxt="n",xlab="",ylab="",yaxs="i")
axis(side=2,at=seq(-1,10,1),las=2)
Ceff = 1.2*10^9
ECS = 3
sigmaQ = 1.0625*10^8
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 4]+1.4, col="grey")
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 8]+2.6, col="grey")
lines(d$year[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff], d[d$sigmaQ==sigmaQ & d$ECS==ECS & d$Ceff==Ceff, 5]+3.4, col="grey")
text(bquote("C=" * .(Ceff/10^9) * "x" * 10^9 * "; " * lambda * "=" * .(round(3.7/ECS,2)) * "; " * sigma[Q] * "=" * .(sigmaQ/10^8) * "x" * 10^8), x=1850, y=4.7, pos=4)

mtext(side=3, expression(Delta * "T + offset (˚C)"), outer=T, line=-3, cex=.8, at=0.05)
mtext(side=1, "Year", outer=T, line=2, cex=.8, at=0.5)

dev.off()
